int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
vector<lower=0>[N] ctry;
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int<lower=0>[N] ctry;
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int<lower=0> ctry;
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
vector[N] ctry;
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
head(sim_ctyr)
head(sim_ctry)
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int ctry[N];
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
N_indivs <- 1000
N_ctries <- 23
b_ses <- 2
b_gini <- -2
b_ses_gini <- 4
u <- rnorm(N_indivs*C,0,6)
sim_ses <- rnorm(N_indivs*N_ctries)
sim_gini <- rep(rnorm(N_ctries), each=N_indivs)
sim_ctry <- rep(seq(1,N_ctries), each=N_indivs)
sim_intercept <- rep(rnorm(N_ctries, 5, 5), each=N_indivs)
ys <- sim_intercept + b_ses*sim_ses + b_gini*sim_gini + b_ses_gini*sim_ses*sim_gini + u
# Adapted from http://www.sanjogmisra.com/blog/stan/
# Adapted from https://github.com/stan-dev/example-models/blob/master/ARM/Ch.3/kidiq_interaction.stan
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int ctry[N];
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
head(sim_dat)
names(sim_dat)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
sim_ctry
length(sim_ctry)
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int<lower=1,upper=N_ctries> ctry[N];
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
seq(1,N_ctries)
class(sim_ctyr)
class(sim_ctry)
length(is.na(sim_ctry))
length(!is.na(sim_ctry))
length(which(is.na(sim_ctry)))
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int<lower=1,upper=N_ctries>[N] ctry;
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int<lower=1,upper=N_ctries> ctry[N];
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=100> sigma;
real<lower=0,upper=100> sigma_ctry_intercepts;
real<lower=0,upper=100> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(100 * mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(df$y),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=df$sim_ctry,
y=df$y,
ses=df$sim_ses,
gini=df$sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
fit <- stan(fit=fit,model_code=model_reg, data=sim_dat, iter=500, chains=2)
fit2 <- stan(fit=fit, data=sim_dat, iter=500, chains=2)
class(fit)
class(fit2)
sim_dat <- list(N=length(ys),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=sim_ctry,
y=ys,
ses=sim_ses,
gini=sim_gini)
fit2 <- stan(fit=fit, data=sim_dat, iter=500, chains=2)
print(fit2)
traceplot(fit2)
?traceplot
traceplot(fit2,pars=c("beta"))
traceplot(fit2,pars=c("mu_beta"))
traceplot(fit2,pars=c("sigma_beta"))
model_reg = "
data {
int<lower=0> N;
int<lower=0> N_indivs;
int<lower=0> N_ctries;
int<lower=1,upper=N_ctries> ctry[N];
vector[N] y;
vector[N] ses;
vector[N] gini;
}
transformed data {
// interaction
vector[N] ses_x_gini;
ses_x_gini <- ses .* gini;
}
parameters {
vector[N_ctries] ctry_intercepts;
vector[3] beta;
real mu_ctry_intercepts;
real mu_beta;
real<lower=0,upper=10> sigma;
real<lower=0,upper=10> sigma_ctry_intercepts;
real<lower=0,upper=10> sigma_beta;
}
transformed parameters {
vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- ctry_intercepts[ctry[i]] + beta[1] * ses[i] + beta[2] * gini[i] + beta[3] * ses_x_gini[i];
}
model {
mu_ctry_intercepts ~ normal(0, 10);
ctry_intercepts ~ normal(mu_ctry_intercepts, sigma_ctry_intercepts);
mu_beta ~ normal(0, 10);
beta ~ normal(mu_beta, sigma_beta);
y ~ normal(y_hat, sigma);
}
"
sim_dat <- list(N=length(ys),
N_indivs=N_indivs,
N_ctries=N_ctries,
ctry=sim_ctry,
y=ys,
ses=sim_ses,
gini=sim_gini)
fit <- stan(model_code=model_reg, data=sim_dat, iter=500, chains=2)
traceplot(fit2,pars=c("sigma_beta"))
traceplot(fit,pars=c("sigma_beta"))
traceplot(fit,pars=c("beta"))
traceplot(fit,pars=c("mu_beta"))
traceplot(fit,pars=c("beta[3]"))
library(stm)
?plot.estimateEffect
library(foreign)
# Use zoo for datetime handling.
library(zoo)
library(lme4)
# Path: (this works for everyone now, but I hadn't moved the data to the dropbox before)
besw = read.dta("~/Dropbox/heytim/data_work_2014/p200506080910_jtw8.dta")
#besw.orig = read.dta("~/Dropbox/heytim/data_work_2014/p200506080910_jtw8.dta")
names(besw)
setwd("~/Documents/tempo/CBEVComparativeReplication/scripts")
library(foreign)
source("./incomedistribution.r")
source("./ca_gen_quintiles_from_distributions.r")
source("./se_gen_quintiles_from_distributions.r")
